MODULE VUPDZ0STL_MOD

USE PARKIND1  , ONLY : JPIM, JPRB

PRIVATE PZ0WN
PUBLIC  VUPDZ0STL

CONTAINS
SUBROUTINE VUPDZ0STL (KIDIA, KFDIA, KLON, KTILES, KSTEP, &
 & KTVL, KTVH, PCVL  , PCVH, PCUR,&
 & PUMLEV5 , PVMLEV5 , &
 & PTMLEV5 , PQMLEV5 , PAPHMS5 , PGEOMLEV5, &
 & PDSN5, &
 & PUSTRTI5, PVSTRTI5, PAHFSTI5, PEVAPTI5 , &
 & PTSKTI5 , PCHAR   , PFRTI   , &
 & PSSDP2  , YDCST   , YDEXC   , YDVEG   , YDFLAKE  , YDURB   ,&
 & PZ0MTI5 , PZ0HTI5 , PZ0QTI5 , PBUOMTI5 , PZDLTI5 , PRAQTI5, &
 & PUMLEV  , PVMLEV  , &
 & PTMLEV  , PQMLEV  , PAPHMS  , PGEOMLEV , &
 & PUSTRTI , PVSTRTI , PAHFSTI , PEVAPTI  , &
 & PTSKTI  , &
 & PZ0MTI  , PZ0HTI  , PZ0QTI  , PBUOMTI  , PZDLTI  , PRAQTI, &
 & PUCURR5 , PVCURR5)  

USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_EXCS  , ONLY : RCHBCD, DRITBL, RCHBBCD, RCHBB, RCHB23A, RITBL, &
 & RCHBA, RCHBHDL, RCDHALF, RCHETB, RCHBD, RCHETA, RCDHPI2, JPRITBL
USE YOS_CST   , ONLY : TCST
USE YOS_EXC   , ONLY : TEXC
USE YOS_VEG   , ONLY : TVEG
USE YOS_FLAKE , ONLY : TFLAKE
USE YOS_URB   , ONLY : TURB
USE YOMSURF_SSDP_MOD

! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     ------------------------------------------------------------------

!**   *VUPDZ0STL*- COMPUTES Z0M,Z0H,Z0Q OVER SEA; SETS Z0H,Z0Q OVER LAND
!                   (Tangent linear)

!     P. Viterbo   ECMWF    22/06/2005     Externalise surf

!     based on

!     J.F. MAHFOUF          E.C.M.W.F.    02/10/95

!     adapted from

!     A.C.M. BELJAARS       E.C.M.W.F.    26/03/90.

!  MODIFIED:
!    M. Janiskova  ECMWF   16/02/2006  modified for tiling of land surface
!                                      (as used in non-linear VUPDZ0 routine)
!    M. Janiskova          21/05/2007  clean-up of z0 initialization
!                                      (as done in vupdz0)
!    G. Balsamo            16/10/2008  lake tile
!    M. Janiskova          21/01/2015  stability param.consistent with NL
!    J. Bidlot             15/12/2018  to use PZ0WN to initialise Z0M over the oceans for step 0
!    J. Bidlot             15/02/2021  Sea state effect in Z0H and Z0Q over the oceans
!    J. McNorton           24/08/2022  urban tile
!    I. Ayan-Miguez (BSC)  Sep 2023    Added PSSDP2 object for surface spatially distributed parameters
!    P. Lopez              July 2025   Added ocean currents (trajectory only)


!     PURPOSE
!     -------

!     DERIVE Z0M,Z0H AND Z0Q FROM SURFACE FLUXES OVER SEA, SET Z0H AND
!     Z0Q OVER LAND AND DERIVE THE BUOYANCY FLUX.
!     (THE T-1 VALUES ARE UPDATED WITH FLUXES FROM THE PREVIOUS TIME
!      STEP)

!     INTERFACE
!     ---------

!     *VUPDZ0STL* IS CALLED BY *SURFEXCDRIVERSTL*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START OF LOOPS
!     *KFDIA*        END OF LOOPS
!     *KLON*         NUMBER OF POINTS IN PACKET
!     *KTILES*       NUMBER OF TILES
!     *KSTEP*        TIME STEP INDEX
!     *KTVL*         LOW VEGETATION TYPE
!     *KTVH*         HIGH VEGETATION TYPE

!     INPUT PARAMETERS (REAL):

!  Trajectory  Perturbation  Description                               Unit
!  PCVL        ---           LOW VEGETATION COVER (CLIMATOLOGICAL)
!  PCVH        ---           HIGH VEGETATION COVER (CLIMATOLOGICAL)
!  PCUR        ---           URBAN COVER
!  PUMLEV5     PUMLEV        WIND X-COMPONENT AT T-1                   m/s
!  PVMLEV5     PVMLEV        WIND Y-COMPONENT AT T-1                   m/s
!  PTMLEV5     PTMLEV        TEMPERATURE AT T-1                        K
!  PQMLEV5     PQMLEV        SPECIFIC HUMUDITY AT T-1                  kg/kg
!  PAPHMS5     PAPHMS        PRESSURE AT T-1                           Pa
!  PDSN5       ---          Total snow depth (m) 
!  PUCURR5     ---           Ocean current U-component (optional)      m/s
!  PVCURR5     ---           Ocean current V-component (optional)      m/s
!  PGEOMLEV5   PGEOMLEV      GEOPOTENTIAL T-1                          m2/s2
!  PUSTRTI5    PUSTRTI       X-STRESS                                  N/m2
!  PVSTRTI5    PVSTRTI       Y-STRESS                                  N/m2
!  PAHFSTI5    PAHFSTI       SENSIBLE HEAT FLUX                        W/m2
!  PEVAPTI5    PEVAPTI       MOISTURE FLUX                             kg/m2/s
!  PTSKTI5     PTSKTI        SURFACE TEMPERATURE                       K
!  PCHAR       ---           "EQUIVALENT" CHARNOCK PARAMETER           -
!  PFRTI       ---           TILE FRACTION                             -



!     OUTPUT PARAMETERS (REAL):

!  Trajectory  Perturbation  Description                               Unit
!  PZ0MTI5     PZ0MTI        NEW AERODYNAMIC ROUGHNESS LENGTH          m
!  PZ0HTI5     PZ0HTI        NEW ROUGHNESS LENGTH FOR HEAT             m
!  PZ0QTI5     PZ0QTI        NEW ROUGHNESS LENGTH FOR MOISTURE         m
!  PBUOMTI5    PBUOMTI       BUOYANCY FLUX                             ????
!  PZDLTI5     PZDLTI        Z/L AT LOWEST MODEL LEVEL 
!  PRAQTI5     PRAQTI        PRELIMINARY AERODYNAMIC RESISTANCE        
!                            FOR MOISTURE

!     METHOD
!     ------

!     SEE DOCUMENTATION

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTILES
INTEGER(KIND=JPIM),INTENT(IN)    :: KSTEP
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVL(KLON)
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVH(KLON)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCVL(KLON)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCVH(KLON)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCUR(KLON)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHMS5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PGEOMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PDSN5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PUSTRTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PVSTRTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PAHFSTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PEVAPTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSKTI5(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCHAR(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PFRTI(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSDP2(:,:)
TYPE(TCST)        ,INTENT(IN)    :: YDCST
TYPE(TEXC)        ,INTENT(IN)    :: YDEXC
TYPE(TVEG)        ,INTENT(IN)    :: YDVEG
TYPE(TFLAKE)      ,INTENT(IN)    :: YDFLAKE
TYPE(TURB)        ,INTENT(IN)    :: YDURB
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PZ0MTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PZ0HTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PZ0QTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBUOMTI5(:,:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PZDLTI5(:,:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRAQTI5(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PUMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PVMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHMS(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PGEOMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PUSTRTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PVSTRTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PAHFSTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PEVAPTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSKTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PZ0MTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PZ0HTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PZ0QTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBUOMTI(:,:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PZDLTI(:,:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PRAQTI(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)  ,OPTIONAL :: PUCURR5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)  ,OPTIONAL :: PVCURR5(:)

INTEGER(KIND=JPIM) :: JL, JTILE

REAL(KIND=JPRB) :: ZUST(KLON,KTILES), ZUST2(KLON,KTILES)
REAL(KIND=JPRB) :: ZRHO(KLON), ZDU2(KLON), ZDUA(KLON)
REAL(KIND=JPRB) :: ZUST5(KLON,KTILES), ZUST25(KLON,KTILES)
REAL(KIND=JPRB) :: ZRHO5(KLON), ZDU25(KLON), ZDUA5(KLON), ZDIV65(KLON)

REAL(KIND=JPRB) :: ZAS, ZAS5, ZCON2, ZIPBL,&
 & ZROWQ, ZROWQ5, ZROWT, ZROWT5, ZUST4, ZUST45, ZWST2, ZWST25,&
 & ZZCDN, ZZCDN5
REAL(KIND=JPRB) :: ZNLEV, ZXLNQ, Z1DZ0Q, ZPRH1, ZPRQ0, ZPRQ
REAL(KIND=JPRB) :: Z0S, Z1S, Z2S, Z3S, Z4S, Z5S, Z6S, Z7S, Z8S
REAL(KIND=JPRB) :: ZEXP1, ZEXP2, ZLOG1
REAL(KIND=JPRB) :: ZNLEV5, ZXLNQ5, Z1DZ0Q5, ZPRH15, ZPRQ05, ZPRQ5
REAL(KIND=JPRB) :: Z0S5, Z1S5, Z2S5, Z3S5, Z4S5, Z5S5, Z6S5, Z7S5, Z8S5
REAL(KIND=JPRB) :: ZEXP15, ZEXP25, ZLOG15
REAL(KIND=JPRB) :: ZDIV15, ZDIV25, ZDIV45, ZDIV55
REAL(KIND=JPRB) :: ZDIV75, ZDIV85, ZDIV95
REAL(KIND=JPRB) :: ZRCPDI, ZRGI
REAL(KIND=JPRB) :: Z0H, Z0Q
REAL(KIND=JPRB) :: ZURBF, ZPCVL, ZPCVH, ZPCVB
! snow
REAL(KIND=JPRB) :: ZSNWGHT5(KLON), ZMLOW5, ZHLOW5, ZTMP5

LOGICAL :: LLCURR, LLINIT

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!             INCLUDE ROUGNESS LENGTH FUNCTIONS
!             ------- -------- ------ ---------
#include "fcz0.h"

!     ------------------------------------------------------------------

!*       1.     INITIALIZE CONSTANTS
!               ---------- ----------

IF (LHOOK) CALL DR_HOOK('VUPDZ0STL_MOD:VUPDZ0STL',0,ZHOOK_HANDLE)
ASSOCIATE(RCPD=>YDCST%RCPD, RD=>YDCST%RD, RETV=>YDCST%RETV, RG=>YDCST%RG, &
 & REPDU2=>YDEXC%REPDU2, REPUST=>YDEXC%REPUST, RKAP=>YDEXC%RKAP, &
 & RNUH=>YDEXC%RNUH, RNUM=>YDEXC%RNUM, RNUQ=>YDEXC%RNUQ, RPARZI=>YDEXC%RPARZI, &
 & RZ0ICE=>YDEXC%RZ0ICE, &
 & LEFLAKE=>YDFLAKE%LEFLAKE, LEURBAN=>YDURB%LEURBAN, &
 & RURBZTM=>YDURB%RURBZTM,RURBZTH=>YDURB%RURBZTH, &
 & RVZ0ML2D=>PSSDP2(:,SSDP2D_ID%NRVZ0ML2D), RVZ0MH2D=>PSSDP2(:,SSDP2D_ID%NRVZ0MH2D),&
 & RVZ0HL2D=>PSSDP2(:,SSDP2D_ID%NRVZ0HL2D), RVZ0HH2D=>PSSDP2(:,SSDP2D_ID%NRVZ0HH2D),&
 & RVZ0M_BARE=>YDVEG%RVZ0M_BARE, RVZ0M_SNOW=>YDVEG%RVZ0M_SNOW, &
 & RVZ0H_BARE=>YDVEG%RVZ0H_BARE, RVZ0H_SNOW=>YDVEG%RVZ0H_SNOW, &
 & RCDFC=>YDEXC%RCDFC)

IF (KTILES.LT.8) THEN
  STOP "Wrong number of tiles in VDFUPDZ0STL"
ENDIF

ZCON2 = 2.0_JPRB/3._JPRB
ZRCPDI = 1.0_JPRB/RCPD
ZRGI = 1.0_JPRB/RG

!     PBL HEIGHT FOR W* - EFFECT

ZIPBL = RPARZI

IF (PRESENT(PUCURR5) .AND. PRESENT(PVCURR5)) THEN
  LLCURR = .TRUE.
ELSE
  LLCURR = .FALSE.
ENDIF

LLINIT= ( KSTEP == 0)

!     ------------------------------------------------------------------

!*         2.      STABILITY PARAMETERS AND UPDATING OF SEA SURFACE
!                  ROUGHNESS LENGTHS

DO JL = KIDIA, KFDIA
  ZDIV15 = 1.0_JPRB/(RD*PTMLEV5(JL)*(1.0_JPRB+RETV*PQMLEV5(JL)))
  ZRHO(JL) = PAPHMS (JL)*ZDIV15 &
   & - RD*PAPHMS5(JL)*ZDIV15*ZDIV15 &
   & * ((1.0_JPRB+RETV*PQMLEV5(JL))*PTMLEV(JL)+RETV*PTMLEV5(JL)*PQMLEV(JL))
  ZRHO5(JL) = PAPHMS5(JL)*ZDIV15
  IF (LLCURR) THEN
    Z0S  = 2.0_JPRB * (PUMLEV (JL) * (PUMLEV5(JL) - PUCURR5(JL)) &
      &  +             PVMLEV (JL) * (PVMLEV5(JL) - PVCURR5(JL)))
    Z0S5 = (PUMLEV5(JL) - PUCURR5(JL))**2 + (PVMLEV5(JL) - PVCURR5(JL))**2
  ELSE
    Z0S  = 2.0_JPRB * (PUMLEV (JL) * PUMLEV5(JL) &
       & +             PVMLEV (JL) * PVMLEV5(JL))  
    Z0S5 = PUMLEV5(JL)**2 + PVMLEV5(JL)**2
  ENDIF
  IF (REPDU2 >= Z0S5) THEN
    ZDU2 (JL) = 0.0_JPRB
    ZDU25(JL) = REPDU2
  ELSE
    ZDU2 (JL) = Z0S
    ZDU25(JL) = Z0S5
  ENDIF
  IF (PDSN5(JL)>0.0_JPRB)THEN
    ZTMP5=PDSN5(JL)/0.25_JPRB
    IF (ZTMP5 < 1.0_JPRB)THEN
      ZSNWGHT5(JL)=ZTMP5
    ELSE
      ZSNWGHT5(JL)=1.0_JPRB
    ENDIF
  ELSE
    ZSNWGHT5(JL)=0_JPRB
  ENDIF
ENDDO

!*         3.   ESTIMATE SURF.FL. FOR STEP 0
!*              (ASSUME NEUTRAL STRATIFICATION)

IF (LLINIT) THEN
  DO JL=KIDIA,KFDIA
    ZDUA (JL) = ZDU2(JL)/(2.0_JPRB*SQRT(ZDU25(JL)))
    ZDUA5(JL) = SQRT(ZDU25(JL))
 
!       - Preliminary value for water
    PZ0MTI (JL,1) = 0.0_JPRB
    PZ0MTI5(JL,1) = PZ0WN(ZDUA5(JL), PGEOMLEV5(JL), PCHAR(JL), RG, RNUM, RKAP, &
         & YDEXC%RITZ0WN, YDEXC%RACDZ0WN, YDEXC%RBCDZ0WN, YDEXC%RXEPSZ0WN, &
         & YDEXC%RUSTMINZ0WN, YDEXC%RPCHARMAXZ0WN, YDEXC%RZ0FGZ0WN)
!       - Sea ice
    PZ0MTI (JL,2) = 0.0_JPRB
    PZ0MTI5(JL,2) = PZ0ICE(RZ0ICE,PFRTI(JL,2))
!       - Wet skin
    PZ0MTI (JL,3) = 0.0_JPRB
    PZ0MTI5(JL,3) = PCVL(JL)*RVZ0ML2D(JL)+PCVH(JL)*RVZ0MH2D(JL)&
     & +(1.-PCVL(JL)-PCVH(JL))*RVZ0M_BARE  
!       - Low Vegetation
    PZ0MTI (JL,4) = 0.0_JPRB
    PZ0MTI5(JL,4) = RVZ0ML2D(JL)
!       - Exposed snow
    ZMLOW5=PCVL(JL)/(1.0_JPRB-PCVH(JL))*RVZ0ML2D(JL)+(1.0_JPRB-PCVL(JL)-PCVH(JL))/(1.0_JPRB-PCVH(JL))*RVZ0M_BARE
    PZ0MTI (JL,5) = 0.0_JPRB
    PZ0MTI5(JL,5)=ZSNWGHT5(JL)*RVZ0M_SNOW+(1.0_JPRB-ZSNWGHT5(JL))*ZMLOW5
!       - High vegetation
    PZ0MTI (JL,6) = 0.0_JPRB
    PZ0MTI5(JL,6) = RVZ0MH2D(JL)
!       - Sheltered snow
    PZ0MTI (JL,7) = 0.0_JPRB
    PZ0MTI5(JL,7) = RVZ0MH2D(JL)
!       - Bare soil
    PZ0MTI (JL,8) = 0.0_JPRB
    PZ0MTI5(JL,8) = RVZ0M_BARE
    IF (LEFLAKE) THEN
!       - Preliminary value for lake
      PZ0MTI (JL,9) = 0.0_JPRB
      PZ0MTI5(JL,9) = 1.E-4_JPRB
    ENDIF
    IF (LEURBAN) THEN
      ZURBF = 0.0_JPRB
      ZPCVL = 0.0_JPRB
      ZPCVH = 0.0_JPRB
      ZPCVB = 0.0_JPRB
      ZURBF=PCUR(JL)
      ZPCVL=MAX(PCVL(JL)*(1.0_JPRB-ZURBF),0.0_JPRB)
      ZPCVH=MAX(PCVH(JL)*(1.0_JPRB-ZURBF),0.0_JPRB)
      ZPCVB=MAX(1.0_JPRB-ZPCVL-ZPCVH-PCUR(JL),0.0_JPRB)
      IF (ZURBF.GT.0.0_JPRB) THEN
       PZ0MTI5(JL,3)=ZPCVL*RVZ0ML2D(JL)+ZPCVH*RVZ0MH2D(JL)&
         & +ZPCVB*RVZ0M_BARE+ZURBF*RURBZTM
       PZ0MTI5(JL,5)=(ZPCVL/(1.0_JPRB-ZPCVH))*RVZ0M_SNOW&
         & +(ZPCVB/(1.0_JPRB-ZPCVH))*RVZ0M_SNOW&
         & +(ZURBF/(1.0_JPRB-ZPCVH))*RURBZTM
      ENDIF
      PZ0MTI (JL,10)=0.0_JPRB
      PZ0MTI5(JL,10)=RURBZTM
    ENDIF

  ENDDO

  DO JTILE=1,KTILES
    DO JL=KIDIA,KFDIA
      ZAS  = (PGEOMLEV(JL)*PZ0MTI5(JL,JTILE)-PGEOMLEV5(JL)*PZ0MTI(JL,JTILE))&
       &     /(RG*PZ0MTI5(JL,JTILE)**2)
      ZAS5 = 1.0_JPRB+PGEOMLEV5(JL)/(RG*PZ0MTI5(JL,JTILE))
      ZLOG1  = ZAS/ZAS5
      ZLOG15 = LOG(ZAS5)
      ZZCDN  = -2.0_JPRB*RKAP**2*ZLOG1/ZLOG15**3
      ZZCDN5 = (RKAP/ZLOG15)**2
      PUSTRTI (JL,JTILE)=ZRHO(JL)*PUMLEV5(JL)*ZDUA5(JL)*ZZCDN5 +&
       & ZRHO5(JL)*PUMLEV (JL)*ZDUA5(JL)*ZZCDN5 +&
       & ZRHO5(JL)*PUMLEV5(JL)*ZDUA(JL) *ZZCDN5 +&
       & ZRHO5(JL)*PUMLEV5(JL)*ZDUA5(JL)*ZZCDN
      PUSTRTI5(JL,JTILE)=ZRHO5(JL)*PUMLEV5(JL)*ZDUA5(JL)*ZZCDN5
      PVSTRTI (JL,JTILE)=ZRHO(JL)*PVMLEV5(JL)*ZDUA5(JL)*ZZCDN5 +&
       & ZRHO5(JL)*PVMLEV (JL)*ZDUA5(JL)*ZZCDN5 +&
       & ZRHO5(JL)*PVMLEV5(JL)*ZDUA(JL) *ZZCDN5 +&
       & ZRHO5(JL)*PVMLEV5(JL)*ZDUA5(JL)*ZZCDN
      PVSTRTI5(JL,JTILE)=ZRHO5(JL)*PVMLEV5(JL)*ZDUA5(JL)*ZZCDN5
      PAHFSTI (JL,JTILE)=0.0_JPRB
      PAHFSTI5(JL,JTILE)=0.0_JPRB
      PEVAPTI (JL,JTILE)=0.0_JPRB
      PEVAPTI5(JL,JTILE)=0.0_JPRB
    ENDDO
  ENDDO
ENDIF

!*         4.      STABILITY PARAMETERS AND FREE CONVECTION
!                  VELOCITY SCALE

DO JTILE=1,KTILES
  DO JL=KIDIA,KFDIA
    ZROWQ  = PEVAPTI (JL,JTILE)
    ZROWQ5 = PEVAPTI5(JL,JTILE)
    ZROWT  = PAHFSTI (JL,JTILE)*ZRCPDI
    ZROWT5 = PAHFSTI5(JL,JTILE)*ZRCPDI
    ZDIV25 = 1.0_JPRB/(PTSKTI5(JL,JTILE)*ZRHO5(JL))
    PBUOMTI(JL,JTILE) = -RG*ZDIV25*(RETV*PTSKTI5(JL,JTILE)*ZROWQ+ZROWT) &
     & + RG*ZDIV25*ZDIV25 &
     & * (PTSKTI5(JL,JTILE)*(RETV*ZROWQ5*PTSKTI5(JL,JTILE)+ZROWT5)*ZRHO(JL) &
     & + ZROWT5*ZRHO5(JL)*PTSKTI(JL,JTILE))
    PBUOMTI5(JL,JTILE) = -RG*(RETV*ZROWQ5*PTSKTI5(JL,JTILE)+ZROWT5)*ZDIV25
    ZUST4  = 2*PUSTRTI(JL,JTILE)*PUSTRTI5(JL,JTILE) &
     & + 2*PVSTRTI(JL,JTILE)*PVSTRTI5(JL,JTILE)
    ZUST45 = PUSTRTI5(JL,JTILE)**2+PVSTRTI5(JL,JTILE)**2
    Z1S5 = SQRT(ZUST45)
    Z1S  =  0.5_JPRB*ZUST4/Z1S5
    ZDIV45 = 1.0_JPRB/ZRHO5(JL)
    ZUST2 (JL,JTILE) = (Z1S *ZRHO5(JL)-Z1S5*ZRHO(JL))*ZDIV45*ZDIV45
    ZUST25(JL,JTILE) =  Z1S5*ZDIV45

!     APPLY W* CORRECTION

    IF (PBUOMTI5(JL,JTILE) > 0.0_JPRB) THEN
      ZWST2  = ZCON2*ZIPBL*PBUOMTI(JL,JTILE) &
       & *(ZIPBL*PBUOMTI5(JL,JTILE))**(ZCON2-1.0_JPRB)
      ZWST25 = (ZIPBL*PBUOMTI5(JL,JTILE))**ZCON2
      ZUST2 (JL,JTILE) = ZUST2 (JL,JTILE)+RCDFC*ZWST2
      ZUST25(JL,JTILE) = ZUST25(JL,JTILE)+RCDFC*ZWST25
    ENDIF

    Z2S5 = SQRT(ZUST25(JL,JTILE))
    IF (Z2S5 >= REPUST) THEN
      ZUST (JL,JTILE) = 0.5_JPRB*ZUST2(JL,JTILE)/Z2S5
      ZUST5(JL,JTILE) = Z2S5
    ELSE
      ZUST (JL,JTILE) = 0.0_JPRB
      ZUST5(JL,JTILE) = REPUST
    ENDIF
    Z3S  = -RKAP*PGEOMLEV (JL)
    Z3S5 = -RKAP*PGEOMLEV5(JL)
    ZDIV55 = 1.0_JPRB/ZUST5(JL,JTILE)
    PZDLTI(JL,JTILE) = ZRGI*(Z3S*PBUOMTI5(JL,JTILE)+Z3S5*PBUOMTI(JL,JTILE)) &
     & * (ZDIV55**3) &
     & - 3._JPRB*ZRGI*ZUST(JL,JTILE)*(Z3S5*PBUOMTI5(JL,JTILE))*(ZDIV55**4)  
    PZDLTI5(JL,JTILE) = ZRGI*(Z3S5*PBUOMTI5(JL,JTILE))*(ZDIV55**3)
  ENDDO
ENDDO

!*         5.    SETTING OF ROUGHNESS LENGTHS
!                ----------------------------

JTILE=1
!   - Ocean open water
DO JL=KIDIA,KFDIA
  ZDIV65(JL) = 1.0_JPRB/ZUST5(JL,JTILE)
  PZ0MTI (JL,JTILE) = -RNUM*ZDIV65(JL)*ZDIV65(JL)*ZUST(JL,JTILE) &
   & + (ZRGI*PCHAR(JL))*ZUST2(JL,JTILE)
  PZ0MTI5(JL,JTILE) =  RNUM*ZDIV65(JL)+(ZRGI*PCHAR(JL))*ZUST25(JL,JTILE)
  Z0H               = RNUH*ZDIV65(JL)
  PZ0HTI5(JL,JTILE) = SQRT(Z0H*PZ0MTI5(JL,JTILE))
  PZ0HTI (JL,JTILE) = 0.5_JPRB * Z0H * (-ZUST(JL,JTILE)*ZDIV65(JL)*PZ0MTI5(JL,JTILE) &
                  & + PZ0MTI (JL,JTILE)) / PZ0HTI5(JL,JTILE)

  Z0Q               = RNUQ*ZDIV65(JL)
  PZ0QTI5(JL,JTILE) = SQRT(Z0Q*PZ0MTI5(JL,JTILE))
  PZ0QTI (JL,JTILE) = 0.5_JPRB * Z0Q * (-ZUST(JL,JTILE)*ZDIV65(JL)*PZ0MTI5(JL,JTILE) &
                  & + PZ0MTI (JL,JTILE)) / PZ0QTI5(JL,JTILE)
ENDDO

JTILE = 2
!   - Sea ice
DO JL=KIDIA,KFDIA
  PZ0MTI (JL,JTILE) = 0.0_JPRB
  PZ0MTI5(JL,JTILE) = PZ0ICE(RZ0ICE,PFRTI(JL,JTILE))
  PZ0HTI (JL,JTILE) = 0.0_JPRB
  PZ0HTI5(JL,JTILE) = RZ0ICE
  PZ0QTI (JL,JTILE) = 0.0_JPRB
  PZ0QTI5(JL,JTILE) = RZ0ICE
ENDDO

JTILE = 3
!   - Wet skin
DO JL=KIDIA,KFDIA
  PZ0MTI (JL,JTILE) = 0.0_JPRB
  PZ0MTI5(JL,JTILE) = PCVL(JL)*RVZ0ML2D(JL)+PCVH(JL)*RVZ0MH2D(JL) &
   & +(1.-PCVL(JL)-PCVH(JL))*RVZ0M_BARE
  PZ0HTI (JL,JTILE) = 0.0_JPRB
  PZ0HTI5(JL,JTILE) = PCVL(JL)*RVZ0HL2D(JL)+PCVH(JL)*RVZ0HH2D(JL) &
   & +(1.-PCVL(JL)-PCVH(JL))*RVZ0H_BARE
  PZ0QTI (JL,JTILE) = PZ0HTI (JL,JTILE)
  PZ0QTI5(JL,JTILE) = PZ0HTI5(JL,JTILE)
ENDDO

JTILE = 4
!   - Low Vegetation
DO JL=KIDIA,KFDIA
  PZ0MTI (JL,JTILE) = 0.0_JPRB
  PZ0MTI5(JL,JTILE) = RVZ0ML2D(JL)
  PZ0HTI (JL,JTILE) = 0.0_JPRB
  PZ0HTI5(JL,JTILE) = RVZ0HL2D(JL)
  PZ0QTI (JL,JTILE) = PZ0HTI (JL,JTILE)
  PZ0QTI5(JL,JTILE) = PZ0HTI5(JL,JTILE)
ENDDO

JTILE = 5
!   - Exposed snow
DO JL=KIDIA,KFDIA
  ZMLOW5=PCVL(JL)/(1.0_JPRB-PCVH(JL))*RVZ0ML2D(JL)+(1.0_JPRB-PCVL(JL)-PCVH(JL))/(1.0_JPRB-PCVH(JL))*RVZ0M_BARE
  ZHLOW5=PCVL(JL)/(1.0_JPRB-PCVH(JL))*RVZ0HL2D(JL)+(1.0_JPRB-PCVL(JL)-PCVH(JL))/(1.0_JPRB-PCVH(JL))*RVZ0H_BARE
  PZ0MTI(JL,JTILE)=0.0_JPRB
  PZ0MTI5(JL,JTILE)=ZSNWGHT5(JL)*RVZ0M_SNOW+(1.0_JPRB-ZSNWGHT5(JL))*ZMLOW5
  PZ0HTI(JL,JTILE)=0.0_JPRB
  PZ0HTI5(JL,JTILE)=ZSNWGHT5(JL)*RVZ0H_SNOW+(1.0_JPRB-ZSNWGHT5(JL))*ZHLOW5
  PZ0QTI(JL,JTILE)=PZ0HTI(JL,JTILE)
  PZ0QTI5(JL,JTILE)=PZ0HTI5(JL,JTILE)
ENDDO

JTILE = 6
!   - High vegetation
DO JL=KIDIA,KFDIA
  PZ0MTI (JL,JTILE) = 0.0_JPRB
  PZ0MTI5(JL,JTILE) = RVZ0MH2D(JL)
  PZ0HTI (JL,JTILE) = 0.0_JPRB
  PZ0HTI5(JL,JTILE) = RVZ0HH2D(JL)
  PZ0QTI (JL,JTILE) = PZ0HTI (JL,JTILE)
  PZ0QTI5(JL,JTILE) = PZ0HTI5(JL,JTILE)
ENDDO

JTILE = 7
!   - Sheltered snow
DO JL=KIDIA,KFDIA
  PZ0MTI (JL,JTILE) = 0.0_JPRB
  PZ0MTI5(JL,JTILE) = RVZ0MH2D(JL)
  PZ0HTI (JL,JTILE) = 0.0_JPRB
  PZ0HTI5(JL,JTILE) = RVZ0HH2D(JL)
  PZ0QTI (JL,JTILE) = PZ0HTI (JL,JTILE)
  PZ0QTI5(JL,JTILE) = PZ0HTI5(JL,JTILE)
ENDDO

JTILE = 8
!   - Bare soil
DO JL=KIDIA,KFDIA
  PZ0MTI (JL,JTILE) = 0.0_JPRB
  PZ0MTI5(JL,JTILE) = RVZ0M_BARE
  PZ0HTI (JL,JTILE) = 0.0_JPRB
  PZ0HTI5(JL,JTILE) = RVZ0H_BARE
  PZ0QTI (JL,JTILE) = PZ0HTI (JL,JTILE)
  PZ0QTI5(JL,JTILE) = PZ0HTI5(JL,JTILE)
ENDDO

IF (LEFLAKE) THEN
  JTILE = 9
!   - Lake
  DO JL=KIDIA,KFDIA
    PZ0MTI (JL,JTILE) = 0.0_JPRB
    PZ0MTI5(JL,JTILE) = 1.E-4_JPRB
    PZ0HTI (JL,JTILE) = 0.0_JPRB
    PZ0HTI5(JL,JTILE) = 1.E-4_JPRB
    PZ0QTI (JL,JTILE) = 0.0_JPRB
    PZ0QTI5(JL,JTILE) = 1.E-4_JPRB
  ENDDO
ENDIF

IF (LEURBAN) THEN
  JTILE = 10
!   - Urban
  DO JL=KIDIA,KFDIA
    ZURBF = 0.0_JPRB
    ZPCVL = 0.0_JPRB
    ZPCVH = 0.0_JPRB
    ZPCVB = 0.0_JPRB
    ZURBF=PCUR(JL)
    ZPCVL=MAX(PCVL(JL)*(1.0_JPRB-ZURBF),0.0_JPRB)
    ZPCVH=MAX(PCVH(JL)*(1.0_JPRB-ZURBF),0.0_JPRB)
    ZPCVB=MAX(1.0_JPRB-ZPCVL-ZPCVH-PCUR(JL),0.0_JPRB)
    IF (ZURBF.GT.0.0_JPRB) THEN
      PZ0MTI5(JL,3) = PZ0MTI5(JL,3)*(1.0_JPRB-ZURBF)+(RURBZTM*ZURBF)
      PZ0HTI5(JL,3) = PZ0HTI5(JL,3)*(1.0_JPRB-ZURBF)+(RURBZTH*ZURBF)
      PZ0QTI5(JL,3) = PZ0QTI5(JL,3)*(1.0_JPRB-ZURBF)+(RURBZTH*ZURBF)
 
      PZ0MTI5(JL,5) = PZ0MTI5(JL,5)*((ZPCVL+ZPCVB)/(1.0_JPRB-ZPCVH))+RURBZTM*(ZURBF/(1.0_JPRB-ZPCVH))
      PZ0HTI5(JL,5) = PZ0HTI5(JL,5)*((ZPCVL+ZPCVB)/(1.0_JPRB-ZPCVH))+RURBZTH*(ZURBF/(1.0_JPRB-ZPCVH))
      PZ0QTI5(JL,5) = PZ0HTI5(JL,5)
    ENDIF
   PZ0MTI (JL,JTILE) = 0.0_JPRB
   PZ0MTI5(JL,JTILE)= RURBZTM
   PZ0HTI (JL,JTILE) = 0.0_JPRB
   PZ0HTI5(JL,JTILE)= RURBZTH
   PZ0QTI (JL,JTILE) = 0.0_JPRB
   PZ0QTI5(JL,JTILE)= RURBZTH
  ENDDO
ENDIF

!*        7.   COMPUTE PRELIMINARY AERODYNAMIC RESISTANCE FOR COMPUTATION
!*             OF APPARENT SURFACE HUMIDITY IN VDFSURF

DO JTILE=1,KTILES
  DO JL=KIDIA,KFDIA
    ZNLEV  = ZRGI*PGEOMLEV (JL)+PZ0MTI (JL,JTILE)
    ZNLEV5 = ZRGI*PGEOMLEV5(JL)+PZ0MTI5(JL,JTILE)
    ZDIV75 = 1.0_JPRB/PZ0QTI5(JL,JTILE)
    Z1DZ0Q  = ZDIV75*ZNLEV-ZNLEV5*ZDIV75*ZDIV75*PZ0QTI(JL,JTILE)
    Z1DZ0Q5 = ZNLEV5*ZDIV75
    ZDIV85 = 1.0_JPRB/Z1DZ0Q5
    ZXLNQ = Z1DZ0Q*ZDIV85
    ZXLNQ5 = LOG(Z1DZ0Q5)
    Z4S  = ZDIV85*PZDLTI(JL,JTILE)-ZDIV85*ZDIV85*PZDLTI5(JL,JTILE)*Z1DZ0Q
    Z4S5 = PZDLTI5(JL,JTILE)*ZDIV85
  
    IF (PZDLTI5(JL,JTILE)  >  0.0_JPRB) THEN
      ZEXP15 = EXP  (-RCHBD*PZDLTI5(JL,JTILE))
      ZEXP1  = -RCHBD*ZEXP15*PZDLTI(JL,JTILE)
      Z5S5 = SQRT(1.0_JPRB+RCHB23A*PZDLTI5(JL,JTILE))
      Z5S  = 0.5_JPRB*RCHB23A*PZDLTI(JL,JTILE)/Z5S5
      ZPRH1  = -RCHBB*(PZDLTI5(JL,JTILE)-RCHBCD)*ZEXP1 &
       & - (RCHBB*ZEXP15+RCHB23A*Z5S5)*PZDLTI(JL,JTILE) &
       & - (1.0_JPRB+RCHB23A*PZDLTI5(JL,JTILE))*Z5S
      ZPRH15 = -RCHBB*(PZDLTI5(JL,JTILE)-RCHBCD)*ZEXP15 &
       & - (1.0_JPRB+RCHB23A*PZDLTI5(JL,JTILE))*Z5S5-RCHBBCD+1.0_JPRB
      ZEXP25 = EXP  (-RCHBD*Z4S5)
      ZEXP2  = -RCHBD*ZEXP25*Z4S
      Z6S5 = SQRT(1.0_JPRB+RCHB23A*Z4S5)
      Z6S  = 0.5_JPRB*RCHB23A*Z4S/Z6S5
      ZPRQ0  = -RCHBB*(Z4S5-RCHBCD)*ZEXP2 &
       & - (RCHBB*ZEXP25+RCHB23A*Z6S5)*Z4S - (1.0_JPRB+RCHB23A*Z4S5)*Z6S
      ZPRQ05 =  -RCHBB*(Z4S5-RCHBCD)*ZEXP25 &
       & - (1.0_JPRB+RCHB23A*Z4S5)*Z6S5-RCHBBCD+1.0_JPRB
    ELSE
      Z7S5 = SQRT(1.0_JPRB-RCDHALF*PZDLTI5(JL,JTILE))
      Z7S  = -0.5_JPRB*RCDHALF*PZDLTI(JL,JTILE)/Z7S5
      ZPRH1  = 2.0_JPRB*Z7S/(1.0_JPRB+Z7S5)
      ZPRH15 = 2.0_JPRB*LOG((1.0_JPRB+Z7S5)*0.5_JPRB)
      Z8S5 = SQRT(1.0_JPRB-RCDHALF*Z4S5)
      Z8S  = -0.5_JPRB*RCDHALF*Z4S/Z8S5
      ZPRQ0  = 2.0_JPRB*Z8S/(1.0_JPRB+Z8S5)
      ZPRQ05 = 2.0_JPRB*LOG((1.0_JPRB+Z8S5)*0.5_JPRB)
    ENDIF
    ZPRQ  = ZXLNQ -ZPRH1 + ZPRQ0
    ZPRQ5 = ZXLNQ5-ZPRH15+ZPRQ05
    ZDIV95 = 1.0_JPRB/(ZUST5(JL,JTILE)*RKAP)
    PRAQTI (JL,JTILE) = ZDIV95*(ZXLNQ-ZPRH1+ZPRQ0) &
     & - RKAP*ZDIV95*ZDIV95*(ZXLNQ5-ZPRH15+ZPRQ05)*ZUST(JL,JTILE)
    PRAQTI5(JL,JTILE) = (ZXLNQ5-ZPRH15+ZPRQ05)*ZDIV95
  ENDDO
ENDDO
END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('VUPDZ0STL_MOD:VUPDZ0STL',1,ZHOOK_HANDLE)
END SUBROUTINE VUPDZ0STL

! Function to compute Z0M for neutral wind conditions
#include "fcz0wn.h"

END MODULE VUPDZ0STL_MOD
